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Abstract 

Within the lowest-order Born approximation, we present an exact calculation of the time dy- 
namics of the spin-boson model in the ohmic regime. We observe non-Markovian effects at zero 
temperature that scale with the system-bath coupling strength and cause qualitative changes in 
the evolution of coherence at intermediate times of order of the oscillation period. These changes 
could significantly affect the performance of these systems as qubits. In the biased case, we find a 
prompt loss of coherence at these intermediate times, whose decay rate is set by y/a, where a is 
the coupling strength to the environment. We also explore the calculation of the next order Born 
approximation: we show that, at the expense of very large computational complexity, interesting 
physical quantities can be rigorously computed at fourth order using computer algebra, presented 
completely in an accompanying Mathematica file. We compute the 0(a) corrections to the long 
time behavior of the system density matrix; the result is identical to the reduced density matrix of 
the equilibrium state to the same order in a. All these calculations indicate precision experimental 
tests that could confirm or refute the validity of the spin-boson model in a variety of systems. 



I. INTRODUCTION 



Novel solid state devices that can control spin degrees of freedom of individual 
electron^ 3, or discrete q uantu m states in subducting circuits^ fl fl Q, sta w 
promise in realizing the ideal of the completely controllable two-state quantum system, 
weakly coupled to its environment, that is the essential starting point for qubit operation in 
quantum computation. From a fundamental point of view, these experimental successes also 
bring us close to embodying the ideal test of quantum coherence as envisioned by Leggett 
many years ago[]|, in which a simple quantum system is placed in a known initial state, is 
allowed to evolve for a definite time t under the action of its own Hamiltonian and under 
the influence of decoherence from the environment, and is then measured. 

Recent exponents, starting witn Q, show that this idea, test can be indented in 
practice. The decay of quantum oscillations due to environmental decoherence is nowP, U, 
|5|, |y] sufficiently weak that some tens of coherent oscillations can be observed. If quantum 
computation is to become a reality, it is believed .g] that these systems will eventually need 
to achieve even lower levels of decoherence, such that thousands or tens of thousands of 
coherent oscillations could be observed. This prospect of producing experiments with ultra- 
long coherence times in quantum two state systems offers a new challenge for theoretical 
modelling of decoherence. Despite the many years of workJj],[n3] following on Leggett's initial 
proposals, there has never been a full, systematic analysis of the most popular description of 
these systems, the spin-boson model, in the limit of very weak coupling to the environment. 

In this paper we provide an exact analysis of the weak coupling limit of the spin boson 
model for the ohmic heat bath, and in the low temperature limit. In this limit the Born 
approximation (to the self energy) should become essentially exact, and we make no other 
approximations in our solutions — in par ticular, no Markov approximation is made. As 
other workers have recently emphasized 3, Q] , understanding the details of the short-time 
dynamics of this model is especially crucial for the operation of these systems as qubits. 

We find important, new, non-Mar kovian effects in this regime. At lowest order in the 
Born expansion of the self energy superoperator, the time dynamics of the model rigorously 
separates into a sum of strictly exponential pieces (the usual "Ti" and "T 2 " decays of the 
Bloch-Redfield model) plus two distinct non-exponential pieces that arise, technically speak- 
ing, from two different kinds of branch cuts in the Laplace transform of the solution of the 
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generalized master equation that we obtain. 

These two contributions both have power-law forms at long times, t > T\,Ti, and thus 
formally dominate the exponentially-decaying parts. But more interesting is that they both 
give new structure to the time evolution at intermediate times t, l/u c < t < T 1; T 2 ; this 
structure typically occurs for t on the order of the oscillation period. (Here, u c is a high 
frequency cut-off of the bath modes, defining the very short time regime, t < l/u c , which is 
of no interest here.) We can explain our results in the language of the double-well potential, 
where the two quantum states are "left" and "right" (L/R), the t = state is pure L, and 
the system oscillates in time via tunneling from L to R. The first branch-cut contribution 
is most important in the unbiased case (L and R energies degenerate) and it causes the 
system, starting immediately in the first quantum oscillation, to spend more time in the R 
well, that is, the opposite well from the one the system is in initially. The second branch- 
cut contribution, present when the system is biased, adds to the amplitude of the coherent 
oscillation, but dies out after an intermediate time which scales like the inverse square root 
of the interaction strength a with the bath. This prompt loss of coherence, whose amplitude 
is proportional to a, changes qualitatively the picture of the initial decay of coherence that 
is so important for discussions of fault-tolerant quantum computation 13j . 

Finally, we set up the next-order Born approximation and do some initial calculations 
with it. This involves computing the self-energy of the master equation to fourth order in 
the system-bath coupling. At this order the self energy is a sum of thousands of separate 
terms; but we find that it is feasible to compute various quantities of physical interest with 
the aid of Mathematica. As an illustration, we provide a full calculation of the steady 
state system density matrix to order a in the limit of low temperature, which requires the 
fourth-order self energy. Given the enormous complexity of the calculation, we find a very 
simple result for the corrections to steady state; they turn out to be identical to those for the 
thermodynamic equilibrium state calculated to the same order in a. Thus, we are able to 
establish rigorously a very strong form of ergodicity for the spin boson model at this order. 

II. GENERALIZED MASTER EQUATION 

We are interested in studying the time dependence of the system density matrix ps(t) = 
Tr Bp{t) with a time-independent system Hamiltonian, and in the presence of a fixed coupling 
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to an environment. An exact equation for p$ - the generalized master equation (GME) - 
isQ 

p s (t) = -iL sPs (t) - i f dt'Z s (t - t')p s (t'), (1) 

Jo 

£ g (t) = -iTi B L SB e- iQLt L SB p B . (2) 

Here the kernel is the self energy superoperator, the system-bath Hamiltonian is 

written 

H = H S + H SB + H B , (3) 

(S = system, B = bath), the Liouvillian superoperator is defined by L x p = [H x ,p], p B = 
e~P HB /Z, (3 = l/k B T, T is the temperature, and Q is the projection superoperator Q = 
1 — p B Ti B . Eq. is written for the case Tr B Hs B p B = 0, and the total initial state is taken 
to be of the form p(0) = ps(0) <g> p B , for an arbitrary ps(0). Since we are interested in the 
case of weak coupling to the bath, we will consider a systematic expansion in powers of this 
coupling L$ B in the self-energy operator 

Retention of only the lowest order term in this expansion, giving the first Born approxi- 
mation, is obtained 15] by the replacement e~ l ® Lt — > e ~ l Q( L s+LB)t j n Thus, in the 
lowest Born approximation, the self energy becomes 

4 2) (t) = -iTr B LsBe- i ( Ls+LB *L 3 BPB. (4) 

We have used the fact here that the expression is unaffected if the Q superoperator is dropped 
in the exponential. 

We now proceed to solve the GME with no further approximations, [l^j This distinguishes 
our work from previous efforts, in which various other approximations Js^ , 11 
wave, Markov, "non-interacting blips", short time) are made (see, e.g. 
We will find that, in particular, avoidance of the Markov approximation endows the solution 
with qualitatively new features. 

We will work out all our results for the ohmic spin-boson model, for which the Hamiltonian 

is 

A e . . 

H s = —<? x + -or Z} (5) 

H S b = ff^E^l + y, (6) 

n 

H B = J2^nb%. (7) 



(secular, rotating 

0, S H H Q 



Here (T x ,y,z are the Pauli operators; we will use ctq = / (identity operator). Also, b\ and b n 
are the creation and annihilation operators of harmonic oscillator n of the bath. With the 
spectral density defined as 

J(u) = J2c 2 J(u-co n ), (8) 

n 

the "ohmic" case is defined by choosing the coefficients c n and the oscillator frequencies uj n 
such that, in the limit of a continuous spectrum, 

J(u,) = |o;e-/- (9) 

Here u; c is an ultraviolet cutoff frequency. 

The first few steps of the solution of the GME do not depend on the details of this model; 
we need only assume that the system Hilbert space is two dimensional, and the system-bath 
coupling has the bilinear form, Hsb — S <g> X (S (X) is an operator in the system (bath) 
space). Under these general circumstances the GME in the Born approximation can be 
rewritten in an ordinary operator form: 

= -iTr 5 <7 M [#s,ps(t)]- fdt'I^t'), (10) 

J 

W) = WHEWOM* -*'))> ( n ) 
i/=i 

V(t') = Re {C(-t')Tr s a u (-t')[a tM ,S}S(-t')} . (12) 

Here (x) = Tr sxps, and the bath correlation function is 

C(t) = Ti B [XX(t) PB ] = C'(t) + iC"{t). (13) 

C and C" denote the real and imaginary parts of the bath correlator), and, for the spin- 
boson model, X = J2n c n(bh + b n ). The time dependent operators are in the interaction 
picture, i.e., 

E(t) = e ^ Hs+HB)t Ee~^ Hs+HB) \ (14) 

for any operator S. 

The GME in Eq. f!10|) can be written in the matrix form 

(&(t)) = R*(a) + k. (15) 

Here a denotes the vector (a x , a y , a z ) T and convolution is denoted A*B = J ' dt'A(t')B(t—t'). 
When the the system Hamiltonian is chosen as in Eq. ((SI), and the system part of the system- 



bath interaction Hamiltonian is S = a z , then we have 



R(t) 



kit) 



( -fr^t) -e5{t)+%K+(t) ^ 

e8(t) - f Jf+(f) -T v (t) -A5(t) 

A5(t) j 

-ffc-(t), -k~(t), Oj 



with 



E = Ve 2 + A 5 



We have introduced the functions 

TiW = 

r„(*) = 

Kit) = 
k-(t) = 
K(t) = 



4A 2 
4A 2 

4eA 
4A 2 a* 



cos(Et)C'{t), 



E 2 

4eA /•* 



E 2 



sin(Et)C'(t), 

dt'sm{Et')C"{t'), 

f dt'(l-cos(Et'))C"(t' 
Jo 



Eq. (|15|) can be solved in the Laplace domain. Defining the Laplace transform as 



e- st f(t)dt, 



the solutions are, for the "standard" initial conditions (a(t = 0)) = (0,0, zq 



Ms)) 

M s )) 
Ms)) 

N(s) 
D(s) 



N(s) 
DM 

AN{s) z 
s D(s) s 



N{s) E 
— 



D(s) 



k-(s) 



E 
A 



E 
A 



A 



K+(a) k-(s) + i-z + k~(s) ) s + —r 1 (s 



E 2 



A 2 



A 2 



f E 2 



E 
A/ 



(16) 
(17) 
(18) 

(19) 
(20) 
(21) 
(22) 
(23) 

(24) 

(25) 
(26) 
(27) 
(28) 
(29) 



To go further, we need an explicit expression for the bath correlator C (t) . For the spin-boson 
model, the well-known formula is 



Cit) 



duiJ{ui){coth.{(3ui/2) cos{ut) + i sin(u;t)). 



(30) 
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For the ohmic case, Eq. (JUJ), Eq. (jSOjl becomes 

1 — iu> c t 



a 



C(t) = --Re V' 



aur: 



(31) 



/5w c / 2(i + u c t) 2 

where ip' is the derivative of the digamma functional. We are not aware that this simple 
exact formula has appeared previously in the literature. 



III. MARKOVIAN LIMIT 



For discussing the exact solution it is instructive to understand the structure of the 
solution in a Markov approximation. This approximation is obtained by replacing all the 
kernels r 1; T y , K+, k~ , and k~ by their forms near s = 0. For all except k~ , this means 
replacing them by constants; k~ has a 1/s divergence at small s. Then the solutions Eqs. 
(|27j) are rational functions of s. If the poles of these rational functions are located at positions 
Sk in the complex s plane, with residues rk/2iri, then the inverse Laplace transform can be 
written {cr^it)) = J2k r k exp(sji). We indicate here that while the residues do depend on the 
label fj, = x,y, z, thepole positions do not, as is shown by the form of Eqs. (|2T|) . 

As is well known[9], there are four poles at positions 



s 1 = 0, 



S-2 



-r° ± %e. 



(32) 



The first pole describes the long-time asymptote of the solution (stationary state), the second 
the purely exponential, "TY'-type decay (relaxation), and the last two (complex conjugate 
paired) describe an exponentially decaying sinusoidal part, the "T 2 "-type decay of coherent 
oscillations. The expressions for the constants in Eq. (J32j) are, to lowest order in a, given 
by 



andH 



r? = 


(T7rA 2 

Tf 1 = ^-coth(^/2), 


(33) 


r° = 




(34) 


E = E + 5E , 


5E = 5E Lamh + 5E Stark , 


(35) 


££.Lamb __ 


aA2 ( C In^V 
E V E) ' 


(36) 


££.Stark __ 


^(Re ^{iE(5/2-K) - HEp/2n)), 
7 


(37) 



where we have dropped terms of order E/u c and higher, C is the Euler constant, and i/j 
is the digamma function [l^j]. These expressions are straightforwardly derivable, and agree 



with the literature 0], except for the energy shift due to vacuum fluctuations, SE Lamh (which 
contains in general ln(u c /E) and not ln(u; c /A)). 

The residues of these poles are, in the limit a — > 0, 

rf = x oc = -(A/E) tanh(/3£/2), r\ = Voo = 0, rf = Zoo = ~{e/E) tanh(/3£/2), 

A/E 2 - Xoo , rf = 0, rl = e 2 /E 2 - Zoo , (38) 



-eA/2E 2 , r y 3A = -A/2E, rf 4 = A 2 /2E 2 . 



' 3,4 

We note that this Markovian theory satisfies the expected fundamental relation 

/oo 
dt(X(t)X) B (Korringa relation); (39) 
-oo 

also, to lowest order in a, the asymptotic values of (er M (£ — > oo)) go to the Boltzmann 
equilibrium distribution of the system^e.g., = —(e/E) tax\h(/3E/2), unlike in the popular 
"non-interacting blip" approximation 



IV. BRANCH CUTS AT T=0 

We now return to the exact solution, examining it in detail at vanishing temperature 

T = 0. In this case Eq. (|31|) becomes 

qui l-u) 2 c t 2 2 u c t 

CT =° {t) = —(l + u*t*r Cr = o(t)=aWc (l + a; c W ( 4 °) 



and the Laplace transform of C is 



C^ =0 (s) = as/2 (— cos (s) Ci (s) — sin (s) si (§)) 

Cr=o( s ) = — ioi/2 (— uj c + s sin (s) Ci (5) — s cos (s) si (s)) , (41) 



where s = s/u; c 19j. There is an important feature of this correlation function that makes 
the Markov solution q ua,itatively incomplete: while the sine integral si( S ) is an analytic 
function of s, the cosine integral Ci(s) behaves like ln(s) for s — ► 17]. This means that 
C(s) is nonanalytic at s = — it has a branch point there. Thus, the exact solutions (a^s)) 
have extra analytic structure not present in the Markov approximation, and the real-time 
dynamics (a^t)) has qualitatively different features in addition to the pole contributions we 
have just discussed. 
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The s = branch point in C(s) leads the kernels ri(s), Ky(s), and k~(s) to have branch 
points at s — ±iE; the kernels T y (s) and k~(s) have three branch points, at s = and 
s = ±iE. Thus, the solutions to the GME (<J x ,y,z{s)) also have three branch points in the 
complex plane. We find by numerical study that the exact solutions still have four poles 
as before, which, for small a, have nearly (but not exactly) the same pole positions and 
residues as in the Markov approximation. 

Thus, the structure of the solutions in the complex s plane is as shown in Fig. Q^. The 
locations of the branch cuts are chosen for computational convenience, as discussed shortly. 
Given this branch-cut structure, the inverse Laplace transform (the Bromwich integral) is 
evaluated by closing the contour as shown. Thus, the exact inverse Laplace transform can 
be expressed as (t > 0) 

M*)> = ^ / dse st (a,(s)) = -L <f dse st (a,(s)) 

Z711 JC ZTTl JCo 

~lT~ Y< qk dxe qkXt ((a^(q k x + 7] k )) - {a^(q k x -%)))• (42) 

Here q k = e l9k and rj k = r]e^ ek ^ w ^ 2 \ with 77 an infinitesimal positive real number. That is, 
r)k is an infinitesimal displacement perpendicular to the direction of branch cut k. For the 
cut choices we have made, Q\ = 57r/4, 6 2 = ir/2, 6> 3 = 3n/2, p% = 0, and P2 = P3 = E. The 
closed-contour integral in the expression can be written as a sum over the four poles, and 
so gives complex exponential contributions to the solution as in the Markovian case. The 
extra terms, the sum over the three branch cuts, are new and give qualitatively different 
features. The contributions of the second and third branch cuts are complex conjugates of 
each other, so we will be discussing them together. 

The contribution of these cuts to the solution is independent of the detailed positioning 
of the branch cuts, so long as they are not moved across a pole; the choice of the direction 
of bcl is a computational convenience — the apparently most natural choice of this cut 
direction, along the negative real axis, passes it essentially on top of the Ti pole, making 
the evaluation of the branch-cut integral numerically inconvenient. As a check, we find that 
the results we discuss now are indeed independent of the cut direction. 

We will do a detailed study of these branch-cut contributions for (cr z (t)) = z(t). We will 
use the following notation for the branch cut terms in Eq. ()42j) : for "branch cut 1" (bcl), 

1 r°° 

Zbci(t) = — / dxe^ x \(o z {q x x + 771)) - {a z {q x x - 771))), (43) 
Lixi J Pl 
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and for two complex-conjugate cuts denoted together as "branch cut 2" (bc2): 

Zbc2 = -7T- Ik / dxe qkXt ((a z (q k x + rj k )) - (cr z (q k x - %))). (44) 

A. Unbiassed case 

For the unbiased spin-boson case, e = 0, an essentially analytic calculation can be done 
for all contributions; we find that these agree, as expected, with the weak-coupling limit 
of the calculations presented in 0. In this case there is no bc2 contribution, z^diif) = 
for all t. The bcl contribution can be obtained analytically to leading order in a: z(t) = 

Zpoles 

{t)+Z bcl (t), 

z bcl (t) = -a{l - At[Ci(At) sin(At) - si(Ai) cos(At)]}. (45) 

This function, plotted along with the pole contribution in Fig. QJd for the choice of parameters 
shown, has the following features: z& c i (t) is negative for all t, it is monotonically increasing, 
and its long-time behavior is 2& c i(t) ~ — 2a/(At) 2 . Also, Zb c i(t = 0) = —a. 

Let us survey, then, the peculiar features that this branch cut contribution introduces 
into the time response z(t). Visualizing the e = spin-boson model as a symmetric double 
well system coupled to its environment, the bcl piece being negative means that, if the 
system is initially in the left well, it will, in the course of coherently tunnelling back and 
forth, spend more time in the right well! This effect becomes strongest at long time, much 
longer than T2, for in this regime the pole contributions are exponentially small, while the 
bcl contribution decays like a power law. Experimentally it may be hard to see the effect in 
this regime (on account of finite-temperature effects, for example), so it is important to note 
that this memory effect appears already at early times, indicating that already in the first 
couple of coherent oscillations, there will be an excess amplitude in the right-well excursions 
as compared with the left- well excursions, by an amount proportional to a. We judge, on the 
basis of a variety of evidence Q , that the Born approximation should be reliable up to a's of 
order 1 — 2%; thus, experiments that look at coherent oscillations accurately at the percent 
level (which, it seems, will ultimately be necessary for performing quantum computation) 
could readily see this bcl effect. 

We note several other interesting features of our solution for e = 0. Taking into account 
the non-Markovian effects, we can do a more precise calculation of the pole positions and 
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residues (only poles 3 and 4 contribute). We find, for T = 0, T 2 = — Re(s 3 ) = T9,r, where, 
as before T° = air A/2, and the renormalization factor r is given by r = (1 — a)/(/t 2 + 
a 2 vr 2 ) < 1, with k = 1 - 2a(l/2 + C - ln(w c /A)). Further, Im(s 3 ) = £ + SE haxah f, with 
f = (k — a7r 2 /2(C — ln(a; c /A)))/ (k 2 + («7r) 2 ). These expressions are obtained as systematic 
expansions in the small parameters T2/E and 5E/E, and they match a direct numerical 
evaluation of the pole positions very well up to a's of a few percent. For the corresponding 
pole residues we find the simple result in leading order r 3 + r 4 = 1 + a + 0(a 2 ). This would 
be impossible in a Markovian theory, in which z(t = 0) = r 3 + r 4 , so that r 3 + r 4 would 
be exactly 1 to all orders in a. In fact this excess pole residue is exactly what is needed to 
cancel out the initial value of the bcl contribution to z(t). We note that our results for the 
residues differ from the weak-coupling expressions in the literature [^] (we are not aware of 
prior reports on the renormalization factors r and f). 



B. Biassed case 



For the biased model (e 7^ 0) the bc2 contributions become nonzero; we find that they 
give other peculiar non-exponential corrections to the solution z(t), very different from the 
bcl contribution. The previous "NIBA" calculations of []| are inapplicable in this case, and 
our results here are completely new. We can do a nearly analytic evaluation of the bc2 
contribution to Eq. (|42jl: Using Eq. (|2*7j) and expanding to lowest order in a, we find for 
the integrand of the sum of the k = 2 and 3 terms of 



• n 2A 2 b~(u) 

Z bc2 {S = ILJ) W — 2 1 h+( \\2 1 U— I \2 - ( 46 ) 

uj [E l — uj z + b + {uj)Y + b \uy 

Here b(iu ± rf) = b + {uj) ± ib~(u), b(s) = a(d(s) + n(s)(s 2 + E 2 )/A), where d(s) and n(s) 
are given by N(s) = A + an(s) (see Eq. flUD) and D(s) = s 2 + E 2 + ad(s) (see Eq. 
Since b~(u) = for |u;| < E, it is reasonable to expect that b~ will grow linearly as one 
passes onto the branch cut; and, in fact, we find from numerical study that a good ansatz 
is b~(uj) = (E — u))b~(u}), with b~(uj) being a weakly varying, real function of uj/ E. With 
this, for uj of order E, Eq. (jlBj) simplifies to 

A 2 b-(E) 1 
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We find that (g7|) should be valid for u > E + b + (E)/2E. Using (gTj) we can do the branch 
cut integral, which gives (for t < \/{o.Exq) — see Appendix for an alternative approach), 



Here is a constant phase shift that we have not determined explicitly (but see Appendix), 
and the dimensionless constants Xq and X\ are 



Since oc a, these constants are independent of a. The last expression for xq comes from 
an evaluation of b + (E): it is directly related to the energy renormalization in the Markov 
approximation, b + (E) = 2E5E Lamb . 

In Fig. El we show a direct numerical evaluation of z^ii). One can see the decay of the 
oscillatory part, which is logarithmic according to Eq. ([48)1 . Even though the decay is very 
non-exponential, it is reasonable to attempt to characterize this decay by a time scale. Eq. 
(148)1 obviously does not work at t = 0, since it is logarithmically divergent. This is not 
surprising, since our calculation has neglected cutoff effects (dependence on uj c ), so Eq. (|4*Kj) 
is not expected to be correct for t < 1/uj c . However, if we consider "early" time to be the 
first half-period of the coherent oscillation, to = ir/E, then Eq. (j4"H|) should be valid and we 
can use it to characterize the decay by determining the time th at which Zf, c2 (t) decreases to 
half its early-time value, i.e., z^ith) = \zbc2{to)- We obtain 



Surprisingly, th oc \j\fa depends non-analytically on a. This explains the effect that is 
evident in Fig. |2J for small a, t h -C T 2 , that is, on the scale of T 2 , there is a very rapid 
loss of coherence as contributed by bc2. This phenomenon may be called a prompt loss of 
coherence, as it would appear experimentally as a fast initial loss of coherence (from 100% 
to (1 — ca)100%, c being some constant near unity), followed by a much slower, exponential 
decay of coherence on the regular T 2 time scale. 

We make a few final remarks about the bc2 calculation. The absolute size of the bc2 
contribution reaches a maximum near the value of e/A used in Fig. El the relative size of this 
contribution continues to increase as |e|/A increases, so that it eventually becomes much 



Zbc2{t) ~ olx\ log(xoa-Et) cos(Et + </>). 



(48) 



x = \b + (E)\/2aE 2 = \6E\/aE 
x x = A 2 t>-{E)/2aE 3 . 



(49) 
(50) 




(51) 
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larger than the pole contribution (but all contributions to z(t) go to zero as \e\/A — > oo). 
When |e| ~ A, we find that, because of the prompt loss of coherence, there is a deficit in the 
total pole contribution, that is, J2k r k — 1 — 0(a) < 1. Even in the absence of an explicit 
branch cut computation, this deficit signals the prompt loss of coherence, in that it indicates 
that the exponentially decaying contributions to z(t) do no account for all the coherence 
near t = 0. Note that this is opposite to the unbiased case, where, as a result of the bcl 
part, there is an excess pole contribution. 

V. NEXT ORDER CALCULATION: STEADY STATE SOLUTION 

Finally, we present the result of a calculation of the corrections to order a to the steady 
state (long time) solution of the GME. To this order, as we will now show, the spin does not 
go to the Gibbs distribution of the uncoupled system (i.e., at T = 0, the ground-state density 
matrix of the isolated spin). However, the result is consistent with the Gibbs distribution 
of the coupled system, giving good evidence for a strong form of ergodicity, even at T = 0. 

These apparently simple corrections, reported below, require an enormous additional 
calculation, in that they require an evaluation of the next order of the Born series. That 
is, we must take the expansion of the self energy superoperator E to fourth order in L$b- 
The formal expression for £( 4 ) is simple enough to generate: it is well known that the full 
Born series is generated by repeated substitution of the following propagator identity into 
the exact expression for £ (Eq. 3.4.7a of 

e -iQLt = e -i Q L t _ i j dUe-^^-^QLsBe- 1 ^ 011 . (52) 

J 

Here 

L = L s + L B . (53) 
This generates the superoperator expression for £( 4 ): 

Tp\t) = {-if f dh rdtzPTTBLsBe^^^QLsBe^^-^QLsBe^^LsBpB 
Jo Jo 

(54) 

This expression can be simplified with the use of the operator identities 

QL = L , PL SB Q = PL SB . (55) 
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Only one factor of the projection superoperator Q = 1 — P survives: 

S (4) (i) = {-if f dti P ^TrBLsBe-^^^IsBe-^^-^g^e-^HsBpB. (56) 
Jo Jo 

Note also that the projector P has been dropped from the expression; since it is immediately 
followed by a trace over the bath it acts as the identity. We can also write E^ 4 ** in several 
equivalent convenient forms using the identity 

e- iLot L SB = L v{ . t) e- lLot , V{t) = e iLot H SB . (57) 

This gives the following two equivalent forms for the self energy: 

S (4) (t) = {-if j^dh l\t 2 Tr B L v[0) L v{tl _ t) QL v{t2 _ t) L v ^ t)PB e- iLs \ (58) 

^ A \t) = {-if J* dt x £ dt 2 e- lLst Tr B L v{t) L v{tl) QL v(t2) L v(0)PB . (59) 

Equation (J58j) can be used to evaluate corrections to the last term of Eq. (fTUj) : we must add 
to Eq. (|12J1 a term of the form 

ijS® = \^s<r^\t)a v . (60) 

The bath part of these traces require the fourth order bath correlator, which using Wick's 
theorem is, at T = 0, 

2 4 



Tr B [X(t 1 )X(t 2 )X(t 3 )X(t 4 )pB] = x 



1 



(l + 00 c {t 3 - t 2 )f{l + 0O c {t A - h)f (l + LO c {t 3 - tx)) 2 (2 + LO c {U ~ t 2 )f 

1 



+ 



(61) 



{i + u c {t 2 -t 1 )f{i + u c {t i -h))\ 

Eqs. (J58H6H) are the starting point of our next-order calculation of the s = residue, which 
gives the long-time asymptote of the density matrix. Every detail of this calculation is 
presented in the accompanying Mathematica notebook. It can be understood why computer 
algebra is necessary for the completion of this calculation if one considers the complexity 
of the above expressions when written out in ordinary operator form. The four nested 
commutators generated by the Liouvillian produces thousands of distinct terms, which all 
need to be integrated and studied in the limit of uj c / ' E — > oo. 
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To illustrate the complexity of this calculation, we give in Appendix B one example of a 
relatively "simple" intermediate result (the integral form for I^(t)) that is obtained in the 
Mathematica notebook. 

Given the enormity of the calculation, the final result is very simple: 



A 

•^oo T OL 



A 3 / u c \ ( A 3 2A\ 



e eA 2 / , u. 



(62) 

Recall that = exactly in the spin-boson model. In this expression all terms that vanish 
in the limit of uj c /E — > oo have been dropped. Note that as in the SE calculation above, we 
see a mild (logarithmic) divergence with the ultraviolet cutoff; all physical quantities that 
we have calculated at this order have no divergence more severe than this. These results 
differ with the 0(a) hoot results reported in Sec. 21.5.2 of g we ceo offer oo explanation 
for this. There is no obvious way of treating the logarithmic divergences in Xoo and by 
introduction of a renormalized A and e, except in the unbiassed case. Nevertheless, the 
expressions given are perfectly physical (x 2 OQ + z"^ < 1) within the expected limits (lo c » E, 
and a < 1/ln^). 

After obtaining the above results, we separately calculated the equilibrium density matrix, 
i.e., 

Trrx e~^ H 2 8 

in the limit T — * and for large u c . Here Z = Tie~^ H , c x = A, and c z = e. We find that 

lim - In Z = - (E + 5E L * mh + au c ) , (65) 

/3^oo fj 2 

with 5E Lamh from Eq. Then it is a simple calculation to show that the equilibrium 

and steady state solutions actually coincide, i.e., 



(0 r x)e< ? ., ^oc = (<7z)eq.- (66) 



While this result is natural, it should not be considered obvious; it provides a rigorous 
demonstration that, to order a, the system is ergodic in a strong sense. 

We give a few final notes about other quantities that require a calculation of E to the 
E^ 4 ) level. The 0(a) corrections to pole positions, given in an earlier section, are unaffected 
by inclusion of Y,^; however, 0(a) corrections to residues of both pole 1 (s — 0) and pole 
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2 for a x and a z are affected by £( 4 ). 0(a) corrections to a y residues, and a X)Z residues of 
poles 2 and poles 3 and 4 are determined entirely by S^ 2 ^; they do not have contributions 
from S^. 



VI. DISCUSSION 

Naturally, many more regimes could be studied using the present approach. For finite 
temperature the time evolution is very different at long times, but it is essentially the same 
as the T = evolution when t < h/kT. Recently, there has been interest in varying both 



the system 



12] and bath 



221 ] initial conditions, as well as in varying the model of the bath 
density of statesj^. For all these circumstances, the systematic Born expansion procedure 
we report here can be done. It is clear on general grounds that the appearance of branch cut 
contributions will not be restricted to the ohmic model, however, the ohmic case is special 
in that the size branch cut contribution is not governed by any small parameter. For any 
superohmic spectral density of the form J{uj) oc u n at low frequencies (n = 1,2,...), C(t) 
will have a power-law dependence at long time, and thus C(s) will have a branch point at 
s = 0. However, the magnitude of the branch cut contribution in the general case goes like 
l/w™~ . So, non-exponential contributions to the dynamics vanish in the physical limit in 
all these other cases. 

Our hope is that, using the present and further exact calculations of the weak-coupling 
behavior of the spin-boson model, a tool will be made available to permit precision ex- 
periments to test the validity of the model (which, at present, is only phenomenologically 
justified) in various physical situations of present interest in quantum information. A fun- 
damentally correct, experimentally verified theory of the system and its environment should 
ultimately be of great value in finding a satisfactory qubit for the construction of a quantum 
information processor. 
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Note that once the choice S = a z is made for the system part of the system-bath interaction 
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form of Eq. © with some choice of the parameters A and e. 

Note that in this theory there is a very simple relationship between a y and a z : Aa y = & z . 



APPENDIX A: SCALING FORM FOR BRANCH CUT INTEGRALS 

By numerical study we find that the branch cut integrals conform to some simple scaling 
laws for small a. If we write the bcl and bc2 integrals as z bc i{t) = dxe~ qixt Zf, c i(s = q\x) 
and Zb C 2(t) = Re J™ dxe lxt z bc2 (s = ix), then we find that for small a and for s « u c , Zb c i(s) 
can be written in a scaling form 

z bcl (x) = (a/E)f 1 (e/A,x/E). (Al) 



18 



But for bc2 a very different scaling law applies: 



z bc2 (x) = (1/E)f 2 (e/A, (x - E)/aE). 



(A2) 



Here /i j2 are dimensionless, "universal" functions that govern the behavior of the branch cut 
contributions for small a. For bcl, the behavior that the scaling law gives is very simple: Eq. 
()A1|) implies that Zb c i{t) = a/i(e/A, Et), where f% is the Laplace transform of the scaling 
function f\. We might have expected this behavior from Eq. (|45|) . from which we can read 
off the scaling function for e = 0. In fact it appears from numerical studies that f\ hardly 
changes as e is varied, except for an overall scale factor; that is, /i(e/A, Et) « a(e/ A)b(Et). 
We find that the scaling function a(r) > is peaked at r = 0. So, the memory effect 
described above for e = persists for finite e, but becomes smaller. For |e| ~ A the bc2 
contribution, which we will describe now, becomes dominant over the bcl one. 

Returning to Eq. (|A2|) . if we write the Fourier transform of the scaling function as 
— Io° e ZXT f 2 (x)dx and consider its polar form fair) = r 2 (T)e l ^ T \ then we obtain 



This shows that bc2 contributes an oscillatory part to the solution, whose "T 2 " decay is 
determined by the features of the scaling function r 2 . A few more observations about f 2 
(obtained initially from numerical study) reveal some crucial properties of the r 2 function: 
1) ^2(0) = 0; 2) |/2(a;)| has a single maximum at x — xq, where xq is some constant of order 
unity; 3) Most important for the present discussion, for x > x f 2 (x) approaches 1/x, that 
is, f 2 (x) ~ Xi/x, where x\ is another real constant of order unity. Fact 3) implies that, for 
T —> 0, r 2 (r) ~ X\ log(x r). That is, we conclude that at sufficiently short time (actually for 
t < l/(aExo), so a relatively long time), 



as stated in the text. 

APPENDIX B: lj$(t) 

As an example of one of many, many intermediate results worked through in the ac- 
companying Mathematica notebook, we give here the expression for lf$(t) (Eq. (|50|l ). in 



Zbc2(t) = ar 2 (aEt) cos(Et + 2 (a£'t)). 



(A3) 



Zbc2(t) = axi \og(x aEt) cos(Et + 0) 



(A4) 
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"simplified" form: 
Jo Jo 



e cos 



(Et) 



e 2 cos(Et) 



cos(Et) 



cos(Et) 



E 2 [ 


-*-^) 2 (*1"*2-^ 


\ / \ W c / 




e 2 cos(Et) 


e 2 cos(Et) 


E 2 (- 

\ 




E 2 (ti-^) 2 (-t + t 2 -^) 2 




9 / T-i ,\ 

e 2 cos(Et) 


9 / n A 

e 2 cos(Et) 


E 2 (- 

\ 




E 2 (t- ±) 2 (-ti + t 2 - ff 


A z 


/ 7 — 1 j 7 — 1 i 7 — 1 I \ 

cos(Et - Et 1 - Et 2 ) 


A 9 / 7 — 1 j 7 — 1 j 7 — 1 j \ 

A 2 cos(£ t- Eti - Et 2 ) 


2E 2 [ 


f -ti-i-) 2 (t-t2-£ 


2E 2 (ti-^) 2 (t-t 2 -^) 2 


A z 


/ 7 — T t T — 1 i 7 — 1 I \ 

cos(Et - Eti - Et 2 ) 


A 9 / 7 — 1 I 7 — 1 j 7 — 1 j \ 

A 2 cos(£ t — Eti - Et 2 ) 


2E 2 { 


-*-^) 2 (*l-*2-^ 


2 Wt-£) 2 (t 1 -t 2 -£) 2 


A 9 

A 2 


/ T — I i 7 — I i 7 — 1 i \ 

cos(Et - Eti - Et 2 ) 


A 9 / 7 — 1 i 7 — 1 i 7 — 1 i \ 

A 2 cos(£ t- Et x - Et 2 ) 


2E 2 (- 

\ 


-ti-±) 2 (-t + t 2 -f 


2E 2 (t 1 -0(-t + t 2 -0 


A 9 

A 2 


/ 7 — 1 I 7 — 1 i 7 — 1 j \ 

cos(Et - Eti - Et 2 ) 


A 9 / 7 — 1 i 7 — I i 7 — 1 i \ 

A 2 cos(£ t- Eti - Et 2 ) 


2 £ 2 (- 


-*-£) 2 (-*!+*»-£ 


2 ^(*-±) 2 (-*i+**-£) 2 


A' 


cos(Et- Eti + Et 2 ) 


A 2 cos(£t - Eti + Et 2 ) 


2^2 ( 


'-ti-t-fit-tt-f 


2E 2 (ti-^) 2 (t-t 2 -^) 2 


A 5 


cos(Et- Eti + Et 2 ) 


A 2 cos(Et - Eh + Et 2 ) 


2£ 2 ( 


\ ^ w c ) ^2 Wcy 


2E 2 (t-£) 2 (ti-t 2 -£) 2 


A 2 


cos(Et- Et x +Et 2 ) 


A 2 cos(Et - Eti + Et 2 ) 


2 (- 


■*i-i) 2 (-* + ^-£ 


^ 2 {h-0{-t + t 2 -^f 


A 2 


cos(£t - Eti + Et 2 ) 


A 2 cos(Et - Eti + Et 2 ) 


2£ 2 (_ 


^ w c ) ( ^ ^2 aj c/ 


) 2 ^ 2 (t-0(-ti + t 2 -0_ 



(Bl) 



This double integral, and many others, are fully evaluated in the Mathematica notebook, in 
the large u c limit. 
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FIG. 1: (a) Structure of the solutions (a fl (s)} in the complex s plane. The four poles pi, P2, P3, 
and p^, are indicated by crosses; the three branch points at s = 0, ±iE are indicated by solid 
circles, and the three branch cuts chosen, bcl, bc2, and bc3, are indicated by dashed lines. The 
inverse Laplace transform requires an integration along the contour C parallel to the imaginary 
axis. This integral may be evaluated by closing with a contour in the left half plane (Co, the 
Bromwich contour), which lies at infinity except for looping back around each of the branch cuts, 
(b) z po i es (t) and Zb c i(t) for the unbiased case, e = 0, A = 1, u c = 30, T = 0, and a = 0.01. t is in 
units of l/E (i.e., E = 1). 
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FIG. 2: Zp i es {t) and 2:fe C 2(t) for the biased case, illustrating the prompt loss of coherence produced 
by bc2. Here E = 1, e/A = —1.38, u; c = 30, T = 0, and a = 0.01. For these parameters, the time 
scale for the prompt loss of coherence (using Eq. I|51|l) is th = 18.98. th is the time at which the 
envelope of falls to half its value at to = ir/E. This time scale is much shorter than the regular 
exponential decay of coherence in z po i es ; for our parameters, T2 = 204.6. 
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